###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##
###*### ESTIMATE PARAMETERS                             ###*###*###*###*###*###*###*###*##
###*### USING SMM                                       ###*###*###*###*###*###*###*###*##
###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##

#' FROM THE PAPER:
#' Child Care Markets, Parental Labor Supply, and Child Development
#' BY:
#' Samuel Berlinski, Maria Marta Ferreyra, Luca Flabbi and Juan David Martin (2023)

# ___________________________________________________________________________________ ####
# START ####

## External packages
library(lubridate)
library(snowfall)
library(readxl)

rm(list = ls())
gc()

## Directory
folderpath <- '~/EstimationCode/'
setwd(folderpath)

# ___________________________________________________________________________________ ####
# INPUTS ####

## Functions
source('functions.R')
source('simul_market.R')

## Sim data
load('SimData_M0.RData')

## Moments
moments     <- as.data.frame(read_excel('AggregateMoments.xlsx', sheet = 1))
sam_moments <- round(moments$Value, 2)
weights     <- 1/moments$S.E.
weights2    <- 1/moments$S.E.

## Parameters vector structure (initial guess)
#' It does not include the 6 parameters for the fathers' wage offer, 
#' which are estimated outside of the algorithm written in this code.
pnames <- c("ma1", "ma2", "ma3", "sa1", "sa2", "sa3", 
            "ma1s", "ma2s", "sa1s", "sa2s", "mu", "tcc", 
            "g0", "g0s", "gf", "gm", "gmh", "gmc", "gsm", "gL", "gH", "gN", "gR", "re", 
           "mwm", "swm", "mwmh", "swmh", "mwmc", "swmc", "rho", "rhoc", 
           "mcL", "mcH", "scL", "scH")
pp0        <- rep(0, length(pnames))
names(pp0) <- pnames
ratio      <- 1

# ___________________________________________________________________________________ ####
# SETUP ####

## Market given settings
J     <- nrow(data)           # number of hh types
NL    <- 7                    # No. of potential entrants
M     <- 40*NL/J/0.3          # Total market size
Fixed <- c(128.90, 322.25)    # Fixed cost for LQ and HQ, respectively
minw  <- 5.15                 # Minumum wage
Cn    <- minw                 # Opp. cost of non-relative caregiver
Cr    <- 0                    # cost of relative caregiver
I     <- 41                   # daily non-labor income
subs  <- 0.35                 # subsidy discount for single moms
Nrel  <- round(J*0.4)         # No. of hh with relative care available

## Observed heterogeneity
q0  <- data$q0      # mental score age 1
ed  <- data$edu     # mom's edu
edf <- data$edf     # dad's edu
wf  <- data$wf      # dad's wage
sm  <- data$single  # single moms
To  <- 16 + 2*sm    # time endowment

## Dad's wage implied random term
ef  <- log(wf)      
mwf <- mean(ef[sm == 0])
swf <- sd(ef[sm == 0])
ef  <- (ef - mwf)/swf
ef[wf == 0] <- 0

## Random draws
set.seed(123)
ev   <- matrix(rnorm(3*J), J, 3) # preferences
ew   <- matrix(rnorm(J),   J, 1) # wages
ec   <- rowMeans(apply(matrix(rnorm(NL*10000), NL, 10000), 2, sort)) # marginal costs
ec   <- (ec - mean(ec))/sd(ec)
reli <- sample(1:J, size = Nrel, replace = F) # hh with relative care available

## Standarization for wages draws
ew[ed == 1] <- (ew[ed == 1] - mean(ew[ed == 1]))/sd(ew[ed == 1])
ew[ed == 2] <- (ew[ed == 2] - mean(ew[ed == 2]))/sd(ew[ed == 2])
ew[ed == 3] <- (ew[ed == 3] - mean(ew[ed == 3]))/sd(ew[ed == 3])

# ___________________________________________________________________________________ ####
# ESTIMATE PARAMETERS ####

## Estimate by minimizing the simulated moments function, eq. (26)
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 3)
sfExport(list = keep, local = FALSE)
n   <- names(pp0)
ans <- optim(pp0[n], fn = crit_funn, parameters = pp0, 
             nampar  = n, sample_moments = sam_moments, 
             weights = weights, method = 'Nelder-Mead', 
             control = list(reltol = 1e-5, maxit = 500))
sfStop()

## Save the new vector of parameter estimates as the optimal solution to eq. (26)
pp1    <- pp0
pp1[n] <- ans$par
pp1    <- conv(pp1)
pp1    <- constraints(pp1)
pp1    <- conv.inv(pp1)
save(pp1, file = 'Estimates_01.RData')

# ___________________________________________________________________________________ ####
# END ####